Genome-wide association study of agronomic traits related to nitrogen use efficiency in Henan wheat

Background Nitrogen use efficiency (NUE) is closely related to crop yield and nitrogen fertilizer application rate. Although NUE is susceptible to environments, quantitative trait nucleotides (QTNs) for NUE in wheat germplasm populations have been rarely reported in genome-wide associated study. Results In this study, 244 wheat accessions were phenotyped by three NUE-related traits in three environments and genotyped by 203,224 SNPs. All the phenotypes for each trait were used to associate with all the genotypes of these SNP markers for identifying QTNs and QTN-by-environment interactions via 3VmrMLM. Among 279 QTNs and one QTN-by-environment interaction for low nitrogen tolerance, 33 were stably identified, especially, one large QTN (r2 > 10%), qPHR3A.2, was newly identified for plant height ratio in one environment and multi-environment joint analysis. Among 52 genes around qPHR3A.2, four genes (TraesCS3A01G101900, TraesCS3A01G102200, TraesCS3A01G104100, and TraesCS3A01G105400) were found to be differentially expressed in low-nitrogen-tolerant wheat genotypes, while TaCLH2 (TraesCS3A01G101900) was putatively involved in porphyrin metabolism in KEGG enrichment analyses. Conclusions This study identified valuable candidate gene for low-N-tolerant wheat breeding and provides new insights into the genetic basis of low N tolerance in wheat. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-023-09922-0.


Background
Bread wheat (Triticum aestivum L.) serves as a primary source of carbohydrates for over two billion individuals globally [1][2][3][4].Although nitrogen (N) is the most important mineral element in the determination of wheat yield [5], N fertilizer is main N source in wheat growth and development.In the past 50 years, the growth rate of nitrogen application in wheat far exceeds that of wheat yield, especially with a nitrogen fertilizer utilization rate of less than 50% [6].Excessive nitrogen application not only affects ecological environment but also reduces production efficiency [7,8].Therefore, a major challenge in crop breeding today is how to improve production efficiency by increasing nitrogen utilization efficiency (NUE).
To address the above issue, the most effective way is to breed "truly green" wheat varieties with high-NUE.However, plant NUE involves a series of physiological processes, including N sensing, uptake, translocation, assimilation, and remobilization [9], while it is a complex trait composed of numerous interacting genetic and environmental factors.Thus, it is critical to identify as many genes related to NUE as possible in the breeding of NUE wheat varieties.
Recently, genome-wide association studies and linkage analysis have successfully identified a series of quantitative trait nucleotides (QTNs) and quantitative trait loci (QTLs) for NUE-related traits in wheat [10][11][12].However, most previous studies involved small number of SNP markers, small population sizes, and experimental environments deviating from production conditions, resulting in their low values in wheat breeding [13][14][15][16][17][18].Therefore, it is necessary and important to mine NUE genes in wheat germplasms.
Grain yield per unit area (GYA) is a comprehensive reflection of crop productivity and the basis for calculating NUE [19].NUE, defined as grain yield per unit of available N from the soil, can be is composed of nitrogen uptake efficiency (NupE) and nitrogen utilization efficiency (NutE) [20].Effective panicle number per unit area (EPN) is not only an important component of yield, but also one of the most related traits to wheat NUE [21].Plant height (PH) significantly affects crop GYA, and has a significantly genetic correlation with NUE [22,23].At present, NUE related traits during the maturity stage have been widely used in identifying rice NUE related loci [9,24,25].Although Henan Province is the largest wheat producing province in China, the genes for wheat NUE related traits have not been fully mined.
To address the above issues, 244 common wheat accessions in Henan Province were genotyped by 660 K wheat microarray and phenotyped by three NUE-related traits (GYA, PH, and EPN) at maturity under high nitrogen (HN) and low nitrogen (LN) conditions in the field in three successive years.3VmrMLM [26] was used to identify QTNs and QTN-by-environment interactions for the three traits.Around these significant loci, candidate genes were predicted based on RNA-seq analysis.This study provides a basis for the further Investigation to detect causal genes of these traits in the future.

Phenotypic data analysis
Three NUE-related traits, PH, EPN, and GYA, for all the 244 accessions were measured under the HN and LN conditions in the field in 2017-2018, 2018-2019, and 2019-2020 (Table 1 and S1).The results showed that, under the LN treatment, PH, EPN, and GYA were decreased only by 6%, 26%, and 22%, respectively, indicating the effectiveness of the LN treatment (Fig. 1A).To reflect NUE, we focused on the ratio of the trait value under LN to the trait value under HN.PH ratio (PHR), EPN ratio (EPNR), and GYA ratio (GYAR) were then  calculated for each accession in each year.The frequency distributions for all three NUE-related trait ratios showed a continuous distribution across all the three environments and BLUP values, coefficients of variation for the BLUP values in three environments ranged from 17.00% (EPN) to 29.53% (GYA) under the HN treatment and from 19.47% (GYA) to 26.97% (PH) under the LN treatment (Table S2), indicating the suitability of these traits for conducting genome-wide association studies (Fig. 1).

Genotype data, LD decay, and population structure analysis
After quality control screening, a total number of 203,224 polymorphic SNP markers were available for subsequent genome-wide association studies.These SNPs were distributed over all the 21 wheat chromosomes (Fig. 2).The above mentioned high-quality SNP markers were used to calculate linkage disequilibrium (LD) decay.The LD decay was measured as the physical distance at which LD dropped to half the maximum value.The LD decay distance was approximately 2.07 Mb (Fig. S1).
Population genetic structure analysis was based on all the 203,224 high-quality SNP markers for 244 wheat accessions.The results showed that a strong peak value of delta K (ΔK) appeared at K = 2 when the number of subpopulations changed from two to five, suggesting the optimal number of sub-populations as K = 2 (Fig. 3C, B).When K = 2, ancient landraces and released cultivars were separated two subpopulations, 1 (Group I) and 2 (Group II).In group I, 135 were landraces.In group II, 109 were released cultivars (Table S1).These observations were further supported by principal component analysis (PCA) (Fig. 3D) and phylogenetic tree construction (Fig. 3A).The above results indicated that the historic periods of release was the main factor for population structure delineation of wheat germplasms in Henan Province in this study.
Identification of QTNs and QTN-by-environment interactions for three NUE-related traits using 3VmrMLM 3VmrMLM was used to perform genome-wide association studies for all the three NUE-related traits.As a result, 133 QTNs were identified by single environment analysis (Table S3; Figs.4A-D, 5A-D, and 6A-D), while 146 QTNs and one QEI were identified by multi-environment joint analysis (Tables S4-5; Figs.4E, 5E, 6E-F).These loci were found to be significantly associated with NUE-related traits and distributed across all the 21 chromosomes.Among these loci, 33 were stably identified across multiple approaches, multiple environments, and multiple traits, these stable loci were distributed on all the chromosomes except 6A, 1B, 6B, and 5D, their proportions of total phenotypic variation.
explained by each QTN (r 2 ) were from 0.13% to 16.59%, and their -log 10 (P) values ranged from 3.20 to 211.56 (Table 2).In particular, 12 QTNs were identified in previous studies, while the remaining 21 were uncovered for the first time (Table 2).In this study, we focused on these new loci.The QTN AX-110463363 on chromosome 3AS with the largest r 2 , which we designated Qphr3A.2, was repeatedly identified in the PHR in 2018 (r 2 : 14.47%) and the BLUP value (r 2 : 16.59%), indicating the existence Based on the stability of the SNPs with significant associations in each environment and the higher phenotype variation explanations, SNP AX-110463363 with largest effect was selected for subsequent candidate gene prediction.Importantly, genetic improvement simultaneously Fig. 4 The Manhattan plot for grain yield per unit area ratio (GYAR) using software IIIVmrMLM.A-C QTNs detected for GYAR in 2018, 2019 and 2020 using single environment module in software IIIVmrMLM.D QTNs detected for GYAR BLUP values using single environment module in software IIIVmrMLM.E QTNs detected for GYAR across all the three environments using multi-environment joint analysis module in software IIIVmrMLM Fig. 5 The Manhattan plot for effective panicle number ratio (EPNR) using software IIIVmrMLM.A-C QTNs detected for EPNR in 2018, 2019 and 2020 using single environment module in software IIIVmrMLM.D QTNs detected for EPNR BLUP values using single environment module in software IIIVmrMLM.E QTNs detected for EPNR across all the three environments using multi-environment joint analysis module in software IIIVmrMLM Fig. 6 The Manhattan plot for plant height ratio (PHR) using software IIIVmrMLM.A-C QTNs detected for PHR in 2018, 2019 and 2020 using single environment module in software IIIVmrMLM.D QTNs detected for PHR BLUP values using single environment module in software IIIVmrMLM.E QTNs detected for PHR across all the three environments using multi-environment joint analysis module in software IIIVmrMLM.F QTN-by-environment interactions (QEIs) detected for PHR using multi-environment joint analysis module in software IIIVmrMLM   [31] for nitrogen efficiency and plant height in wheat might be carried out by selecting for a single large-effect QTN.

Prediction of NUE-related candidate genes
Based on the LD decay distance of 2.07 Mb in Fig. S1,  2.07 Mb on each side of the significant locus qPHR3A.2for NUE was used to mine NUE-related candidate genes from the Chinese spring version 1.0 genome (http:// 202.194.139.32/).As a result, there were 52 annotated genes in the region (Table S6).Further KEGG enrichment analysis showed that 11 of them could be functionally enriched mainly in the metabolic, genetic information processing and environmental information processing pathways in common wheat, especially, there were seven metabolism-related candidate genes (Fig. 7A and Fig. S2) ( [34][35][36], https:// www.kegg.jp/ kegg/ kegg1.html).
To further determine candidate genes, transcriptome analysis was performed on all the 52 annotated genes using cultivars with different low nitrogen tolerance (Table S7).Results showed that four of them were differentially expressed after induction, i.e., the expression levels were significantly changed in at least one accession with different N treatments, indicating that the four genes were transcriptionally induced by the LN treatment (Fig. 7B; Table S6).Among the four DEGs, two were related in metabolic pathways in the analysis of KEGG and referred as candidates for NUE-related genes (Fig. S2; Table S6).In detail, TraesCS3A01G102200 was differentially expressed only in Jinzihong at 15 DPA (LN-JE vs HN-JE).However, TraesCS3A01G101900 was differentially expressed between the HN and LN treatments both in Bainong 207 and Jinzihong at anthesis (LN-BA vs HN-BA and LN-JA vs HN-JA).Although Fig. 7 Candidate genes around qPHR3A.2.A KEGG enrichment analysis of 52 annotated genes.The pathways enriched for the target genes are in red font.The KEGG image used in this figure is licensed by KEGG copyright and has been changed accordingly.B Fold change of 52 candidate genes.The red, black, and green represent higher in low N than high N of a particular accession at a particular growth stage, close to high N, and lower than high N of a particular gene, respectively.Differentially expressed genes are in red font.Target genes are in blue boxes.LN low N treatment, HN high N treatment, B Bainong 207, J Jinzihong, A flowering stage, E 15 days past anthesis TraesCS3A01G101900 was no differentially expressed in these two varieties at 15 DPA (LN-BE vs HN-BE and LN-JE vs HN-JE), TraesCS3A01G101900 was down-regulated and up-regulated expressed at 15 DPA in Bainong 207 and Jinzihong (LN-BE vs HN-BE and LN-JE vs HN-JE), respectively.Among the two candidates, TraesC-S3A01G101900 was annotated in KEGG as a chlorophyllase and encoded hydrolase superfamily structural domain protein in NCBI conserved structural domain analysis, which was found in recently previous studies to be widely involved in abiotic stress processes [37].In conclusion, TraesCS3A01G101900 was regarded as preferred candidate gene of qPHR3A.2,named TaCLH2.

Discussion
Reliable results and significant progress GWAS for wheat NUE in 244 landraces and released varieties (lines) in Henan Province, China provides new insight into the genetic basis of the important quantitative trait and varietal breeding.This study made reference to the results of Xing et al., 2022, and meet the increasing demand for wheat production, PHR, EPNR, and GYAR at wheat mature stage were used for GWAS.The results of GWAS showed that the reported NUE-related genes/ QTLs (qGNSR2B.3,TaGS1.3,MQTL5B.5, MQTL4D.2) could be detected repeatedly in various environment.Among them, TaGS1.3, which was directly associated with NUE, were repeatedly detected in different environments in this study [12].In summary, this analysis shows that these three traits can be taken as the indicator of NUE.
The 3VmrMLM model achieves very good performance, e.g., high power and accuracy, low false positive rate, and correct detection of all types of QTNs [26,35].RNA-seq has become an effective technique to examine genome-wide gene expression level [38].Nevertheless, a very large set of DEGs is typically obtained, which made it difficult to identify underlying interesting candidate genes.Moreover, combined 3VmrMLM with RNA-seq has been recently used to predict candidate genes in soybean [39], rice [40], maize [41], and others.In the present study, a new comprehensive GWAS method, 3VmrMLM, was used for detecting QTNs, and QEIs of NUE-related traits under various nitrogen treatments.As a result, 12 stable QTNs had been reported in Wheat (Table 2).

Comparison of N-efficient loci with reported loci in mature wheat
To investigate the relationship between the 33 replicated NUE SNP markers detected at maturity and the reported wheat NUE loci, a literature search revealed that 12 of the SNP markers may be closely linked to previously reported NUE loci and 21 markers may represent new NUE loci.Comparison revealed that AX-109315961, which was associated with relative effective spike per unit area and relative seed yield per unit area on chromosome 2AS in this study, might be related to Chr2A_SNP_03717, which has been reported to be associated with seed yield at field maturity.This association occurs because the two loci, SNP marker AX-109315961 and the previously reported Chr2A_SNP_03717, are approximately positioned 479,000 bp from each other on the Chinese spring reference genome [32].SNP AX-109519142 on chromosome 3AS is 574 kb away from the previously reported BobWhite_c18256_105, as a consequence, SNP AX-109519142 is identical to BobWhite_c18256_105 [33].SNP marker AX-109874103 on chromosome 5AL is the same locus as qFsnR5A.2,because AX-109874103 is located between SNP markers AX-110964357 and AX-109942297, which are closely linked to qFsnR5A.2 [27].AX-108917610 on chromosome 5AL is the same locus as the qGNSR5A.2NUE locus, as this locus is located both between the SNPs AX-109948410 and AX-111152617 that are tightly linked to both sides of qGNSR5A.2[27].
AX-111497324 on chromosome 2BL is at the same locus as qGNSR2B.3 because AX-111497324 is located between markers AX-109374416 and AX-109848914 that are linked to both sides of qGNSR2B.3 [27].The SNP AX-109412797 on chromosome 4BL is 83 kb away from BobWhite_c8266_582, and therefore AX-109412797 is considered to be the same locus as BobWhite_c8266_582 [28].SNP AX-110432662 on chromosome 5BS is identical to MQTL5B.5 because SNP AX-110432662 is located between SNP markers AX-110612519 and AX-109473183, which are closely linked to MQTL5B.5 [31].The SNP marker AX-109239751 on chromosome 2DS is 695 kb away from Kukri_c16477_157, and given this, AX-109239751 and Kukri_c16477_157 are also considered to be the same locus [28].SNP AX-89699893 on chromosome 3DL is the same locus as AX-109688668, as SNP AX-89699893 is 1.59 Mb away from AX-109688668 [29].Similarly, SNP AX-110061772 on chromosome 3DL is 0.32 Mb away from the reported SNP AX-109649974, so we consider these 2 loci as the same locus [29].SNP AX-89605902 on chromosome 4DS is tightly linked to the wheat glutamine synthetase gene TaGS1.3,as this locus is only 660 kb away from TaGS1.3 (TraesCS4D02G047400) [30].The SNP marker AX-111684216 on chromosome 4DL is located between SSR Xwmc399, which is closely linked to MQTL4D.2, and SNP AX-109334705, so AX-110242843 is considered to be the same locus as MQTL4D.2[31].

Future prospects
A review of the literature revealed that 12 of the 33 replicate detected QTNs were closely linked to reported QTLs, loci or genes significantly associated with NUE and NUE-related traits in wheat.A novel locus on chromosome 3A, SNP AX-110463363, had an explanation rate between 1.66%-16.59%,indicating that genes affecting NUE are likely to be present on chromosome 3A.Next, we performed gene annotation of the candidate interval qPHR3A.2for the new locus of NUE, and a preferred candidate gene encoding chlorophyllase was detected on chromosome 3A (AX-110463363).The hydrolase, chlorophyllase, is responsible for the breakdown of chlorophyll a, thereby generating dephyllophyll a and phytol which is known to play important roles in plant fruit ripening, leaf senescence, and processes resulting from biotic and abiotic stresses.
It has been suggested that citrus CLH plays a central role in chlorophyll degradation during decolorization of citrus fruits [42].Hu et al. (2015) suggested that CLH is involved in secondary defense of Arabidopsis against herbivorous insects [43].In a study on the effects of nitrogen fertilization on leaf senescence and carbon and nitrogen metabolism in late maize reproduction, Du et al. (2020) found that under suitable nitrogen fertilization conditions, the CLH activity of cob leaves was relatively low and chlorophyll content decreased slowly, which was conducive to the prevention of premature leaf senescence in cob leaves [44].Zhang et al. (2020) found significant up-regulation of expression of genes encoding chlorophyllase in transcriptome analysis of N reduction on leaves of potato varieties with different N efficiency [45].A study by Tian et al. (2021) showed that in Arabidopsis CLH1 plays an important role in young leaf photoprotection by activating chlorophyll degradation activity through N-terminal shear post-translational modification, degrading chlorophyll and assisting chloroplast protease-mediated degradation of PSII core protein [46].

Conclusion
Thirty-three significant QTNs were stably identified across different traits, environments, and approaches using compressed variance component mixed model method, 12 previously reported loci confirmed the reliability of the results in this study, more importantly, 21 novel loci were identified.Around the locus qPHR3A2, TaCLH2 (TracesCS3A01G101900) was found in KEGG and gene expression analysis to be associated with NUE, although its biological function and molecular mechanism required further research in the future.

Plant material
The association panel was comprised of 135 landraces in 1940s and 109 improved cultivars from 1940 to 2010s (main cultivars, backbone parents, and pioneer varieties), and their seeds were provided by Henan Academy of Crops Molecular Breeding.The formal identification of the plant material used in this study was performed by Professor Weigang Xu of Henan Academy of Agricultural Sciences, and the voucher specimen of this material has been deposited in Institute of Crop Molecular Breeding, Henan Academy of Agricultural Sciences, Henan, China.

Field trait measurement and statistical analysis
All the accessions were planted in three growing seasons in the Henan Research and Development Center for Modern Agriculture (Yuanyang County, Henan Province; 35• 00′ N 113• 40′ E, 77 m a.s.l.).All the trials were arranged in a completely randomized block design with three replications.Each plot consisted of four rows with 2 m in length and 0.23 m in width between rows, and sowing was conducted manually at a seeding rate of 130 seeds/m 2 on October 11, 2017, October 10, 2019, and October 10, 2020.Two N treatments were applied: low N treatment (LN) and high N treatment (HN).The basic physicochemical characteristics and fertilization level are detailed in Table 1.
As described in Peng et al. (2021), PH, EPN, and grain yield per unit area (GYA) were measured [47].The method described in Shi et al. (2022) was used to calculate low N tolerance for PH, GYA, EPN, and the best linear unbiased estimation (BLUP) [41].As described in Tang et al. (2019), the processed data post normalization was used as the adjusted phenotypic values for genomewide association studies [9].For the zero mean normalization, the equation was where x Z is the zero mean normalized phenotypic val- ues, x is the relative phenotype values, µ is the phenotype mean, σ is the standard deviation.
Descriptive statistical analysis and one-way ANOVA were performed using IBM SPSS v22.0 (SPSS, Chicago, IL, USA).

Genotyping and marker screening
Total genomic DNAs of 244 accessions were extracted from fresh seedling leaves using a modified CTAB method, following Du et al. (2021) [48].Genotyping for all 244 wheat accessions was carried out using SNP 660 K wheat microarray developed by Zhongyujin Marker Biotechnology Co., Ltd., Beijing, China.To reduce errors, SNP markers with minor allele frequency (MAF) < 0.05, missing data > 0.1 and heterozygosity > 0.1 were removed.

Population structure and linkage disequilibrium (LD) analysis
The population structure of 244 wheat accessions was determined using the STRU CTU RE v2.3.4 software [49].
The number of subgroups (K) was set from 2 to 5. We set BURNIN = 10,000 and NUMREPS = 110,000, respectively, and ran 5 replicates for each value of K.

Genome-wide association studies
Q matrix produced by Structure was included as covariate in the analysis to control for potential populations structure.To control polygenic background, marker inferred kinship matrix K was calculated by the 3VmrMLM software [26].The single environment module of the 3VmrMLM method was used to identify QTNs for three NUE-related traits in each environment, while multi-environment module was used to detect QTNs and QTN-by-environment interactions [50].As described in Li et al. (2022b), the critical P-value and LOD thresholds were set as 0.05/m and 3.0, respectively, for significant and suggestive QTNs, where m is the number of markers [50].

Identification of candidate genes
Based on the LD attenuation distance, the related genes in a flanking region of the physical location of the most significant/suggested loci were annotated using wheat reference genome (CS RefSeq version 1.0) [51].
Bainong 207 (Low-N sensitive) and Jinzihong (lownitrogen-tolerant) were used as the materials of RNAseq analysis, and the field experiments were described as previous ones.At the anthesis stage, the main shoots of the plants was tagged.10 flag leaves of each plot was randomly sampled at anthesis, 15 days past anthesis (DPA).Samples were quick frozen in liquid nitrogen after sampling and stored at − 80 °C.Total RNA extraction, mRNA library preparation, and Illumina sequencing were entrusted to Guangzhou Genedenovo Biotechnology Co., Ltd.(Guangzhou, China).
Differential expression analysis and KEGG enrichment analysis were carried out on with the Omicsmart website (https:// www.omics mart.com) and KEGG databases [34][35][36].Expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) was used to evaluate genes expression level.Only expressed genes with mean FPKM ≥ 1 were considered.False discovery rate (FDR) ≤ 0.05 and |log 2 (fold change)|≥ 1 were considered as screening criteria for differentially expressed genes (DEGs).The heatmap module in TBtools (v1.120) was used to cluster and draw the heat map [52].
All methods were carried out in accordance with relevant guidelines.
• fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year

•
At BMC, research is always in progress.

Learn more biomedcentral.com/submissions
Ready to submit your research Ready to submit your research ?Choose BMC and benefit from: ? Choose BMC and benefit from:

Fig. 2 Fig. 3
Fig. 2 The density distribution of single nucleotide polymorphisms (SNPs) in wheat.Red and gray horizontal bars show genomic regions that are rich and poor in SNPs, respectively

Table 1
[13]soil basal fertility and fertilizer application rate for each treatment in 2017-2018, 2018-2019, and 2019-2020 a Information of the soil basal fertility and fertilizer application rate for each treatment in 2017-2018 and 2018-2019 from Peng et al. (2022)[13]

Year Treatment Determination parameter (kg⋅ha −1 ) Base fertilizer (kg⋅ha −1 ) Topdressing N (kg⋅ha −1 ) Alkali hydrolysable N Available phosphorus Available potassium N P 2 O 5 K 2 O
Phenotypic diversity of 244 wheat accessions.A Frequency distributions and differences for PH, EPN, and GYA of the wheat accessions measured under LN and HN in 2018, 2019, 2020, and the BLUP across all three environments.
PH: plant height, EPN: effective panicle number per unit area, GYA: grain yield per unit area.BLUP: Best linear unbiased predictions.Box edges represent the 0.25 and 0.75 quantiles with the median values shown by bold lines.Whiskers extend by no more than 1.5 times the interquartile range, and the remaining data are indicated by dots.P values were calculated with Student's t test.B Frequency distribution histogram of zero mean normalized phenotypic values of PHR,